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Abstract 

A  newly  developed  hybrid-coordinate  ocean  circulation  model  is  documented  and  tested.  Coordinate 
surfaces  in  this  model  adhere  to  isopycnals  wherever  this  does  not  violate  minimum  layer  thickness  re¬ 
quirements;  elsewhere,  coordinate  surfaces  are  geometrically  constrained.  The  intent  of  this  approach,  some 
of  whose  features  are  reminiscent  of  the  Arbitrary  Lagrangian-Eulerian  (ALE)  technique,  is  to  combine  the 
best  features  of  isopycnic-coordinate  and  fixed-grid  circulation  models  within  a  single  framework.  The 
hybrid  model  is  an  offshoot  of  the  Miami  Isopycnic  Coordinate  Ocean  Model  whose  solutions,  obtained 
under  identical  geographic  and  forcing  conditions,  serve  as  reference.  Century-scale  simulations  on  a 
coarse-mesh  near-global  domain  show  considerable  similarities  in  the  modeled  thermohaline-forced  cir¬ 
culation.  Certain  architectural  details,  such  as  the  choice  of  prognostic  thermodynamic  variables  (p,  S 
versus  T.  S )  and  the  algorithm  for  moving  coordinate  surfaces  toward  their  reference  isopycnals,  are  found 
to  only  have  a  minor  impact  on  the  solution.  Emphasis  in  this  article  is  on  the  numerical  resiliency  of  the 
hybrid  coordinate  approach.  Exploitation  of  the  model’s  flexible  coordinate  layout  in  areas  of  ocean 
physics  where  pure  isopycnic  coordinate  models  only  have  limited  options,  such  as  mixed-layer  turbulence 
parameterization,  will  be  the  subject  of  forthcoming  articles.  ©  2001  Elsevier  Science  Ltd.  All  rights 
reserved. 


1.  Introduction 

Faithful  replication  of  the  conservation  properties  inherent  in  the  laws  of  physics  is  a  near- 
universal  goal  in  numerical  modeling.  In  the  case  of  rotating,  stratified  geophysical  flows,  the 
conservation  properties  of  the  equations  comprised  in  a  particular  model  can  be  manipulated  not 
only  by  changing  the  size  of  the  computational  mesh  or  the  sophistication  of  the  operators 
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mimicking  spatial  and  temporal  differentiation,  but  also  by  transforming  the  conservation  laws 
from  a  Cartesian  to  a  nonCartesian  vertical  coordinate.  The  latter  manipulation  explains  the 
multitude  of  vertical  coordinates  proposed  over  the  years  for  use  in  atmospheric  and  oceanic 
circulation  models. 

Perhaps  the  most  striking  example  of  such  manipulation  of  the  governing  equations  is  the  use 
of  entropy  as  vertical  coordinate  (Montgomery,  1937;  Starr,  1945).  This  transformation  -  made 
possible  in  the  oceanic  context  by  the  fact  that  entropy  (or  its  proxy,  potential  density  p  t)  varies 
monotonically  with  depth  in  most  parts  of  the  world  ocean  -  reduces  the  First  Law  of  Ther¬ 
modynamics  in  the  important  case  of  adiabatic  flow  to  the  simple  statement  ppot  =  dp  Jdt  =  0. 
This  in  turn  renders  adiabatic  fluid  motion  two-dimensional  in  (x,y,p  t)  space.  These  simplifi¬ 
cations  readily  carry  over  to  the  algebraically  approximated  equations  in  numerical  models  using 
entropy  as  vertical  coordinate. 

Fluid  mixing  and  stirring  likewise  are  handled  more  concisely  if  carried  out  in  ppot  space  where 
a  sharp  line  can  be  drawn  between  eddy-driven  isopycnal  stirring  and  diapycnal  mixing.  Both 
processes  are  important  in  the  ocean  but  typically  act  on  vastly  different  space  and  time  scales. 

One  particular  incarnation  of  fluid  laws  formulated  in  (x,y,p  t)  space  is  the  so-called  stacked 
shallow  water  model,  a  system  of  equations  describing  the  dynamics  of  a  stratified  fluid  in  terms  of 
the  dynamics  of  a  set  of  interacting  constant-potential  density  layers.  Stacked  shallow  water 
models  have  been  widely  used  as  theoretical  tools  elucidating  aspects  of  large-scale  oceanic  and 
atmospheric  flow  (example:  baroclinic  instability).  Use  of  this  paradigm  in  numerical  models  also 
has  a  long  history,  at  least  in  process  models  dealing  with  particular  aspects  of  fluid  behavior.  On 
the  other  hand,  full-fledged  oceanic  circulation  models  built  upon  the  use  of  entropy  as  vertical 
coordinate  have  been  developed  only  quite  recently  (Bleck  et  al.,  1992;  Oberhuber,  1993). 

The  tradeoffs  between  advantages  and  disadvantages  of  isopycnal  modeling  (the  word  iso¬ 
pycnal  refers  here  to  the  property  ppot  =  const.)  have  been  the  subject  of  vigorous  debate.  While 
there  is  little  disagreement  about  the  potential  advantages,  the  discomfort  about  the  disadvan¬ 
tages  of  isopycnal  modeling  varies  greatly  among  numerical  practitioners.  The  most  serious 
drawbacks  of  an  entropy-based  vertical  coordinate  are  its  degeneracy  in  unstratified  and  statically 
unstable  water  columns  and,  to  a  lesser  extent,  the  possibility  that  coordinate  surfaces  intersect  the 
sea  surface  or  bottom.  While  robust  numerical  methods  have  been  developed  to  deal  with  the 
latter  aspect,  the  prevailing  opinion  in  the  community  is  that  modeling  of  oceanic  flow  in  un¬ 
stratified  -  let  alone  convectively  unstable  -  regions  cannot  properly  be  done  in  potential  density 
space. 

Another,  more  practical  drawback  of  isopycnal  modeling  is  the  requirement  that  the  vertical 
grid  in  an  isopycnic  model  span  a  wide  enough  range  to  capture  the  extremes  of  water  density 
found  anywhere  in  modeled  domain.  Many  grid  points  in  density  space  are  thereby  “wasted”.  To 
put  it  differently:  The  world  ocean  does  not  fit  snugly  into  a  rectangular  box  in  (x,y,  p)  space. 

This  paper  presents  initial  results  obtained  with  a  numerical  model  which,  while  still  primarily 
an  isopycnic  coordinate  model,  addresses  all  of  the  above  concerns  by  allowing  coordinate  sur¬ 
faces  to  deviate  locally  from  isopycnals  wherever  the  latter  either  fold,  outcrop,  or  generally 
provide  inadequate  vertical  resolution  in  portions  of  the  modeled  domain. 

The  concept  of  “hybrid”  coordinates  that  act  in  the  above-described  manner  is  not  new.  The 
groundwork  for  formulating  the  equations  governing  geophysical  fluid  flow  in  a  generalized 
vertical  coordinate  was  laid  by  Starr  (1945),  Kasahara  (1974),  and  Bleck  (1978a).  Coordinate 
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surfaces  that  adhere  to  isopycnals,  except  in  regions  where  these  outcrop,  were  introduced  two 
decades  ago  by  Bleck  and  Boudra  (1981)  (hereafter  BB81)  in  a  study  of  the  wind-driven  circu¬ 
lation  in  an  idealized  ocean  basin.  The  rationale  for  the  use  of  hybrid  coordinates  at  that  time,  at 
least  from  the  perspective  of  the  BB81  authors,  was  lack  of  numerical  skill  in  dealing  with  the 
outcrop  problem.  Hence,  the  hybrid  coordinate  concept  was  promptly  shelved  when  techniques 
for  numerically  solving  the  shallow-water  equations  in  the  limit  of  zero  layer  thickness  became 
available  (e.g.,  Bleck,  1984;  Bleck  and  Boudra,  1986).  Elsewhere,  however,  hybrid  ocean  modeling 
continued  (Zebiak  and  Cane,  1987;  Schopf  and  Loughe,  1995). 

The  BB81/Schopf  hybridization  approach  differs  from  other  generalized  coordinate  schemes 
developed  over  the  years  (e.g.,  Bleck,  1978b,  and  references  therein;  Zhu  et  al.,  1992;  Konor  and 
Arakawa,  1997)  in  that  it  does  away  with  an  analytic  formula  specifying  where  in  a  fluid  column  a 
given  grid  point  is  located.  Rather,  the  BB81/Schopf  scheme  assigns  each  grid  point  a  reference  or 
“target”  isopycnal  and  continually  tries  to  move  individual  grid  points  that  have  become  sepa¬ 
rated  from  their  target  isopycnal  back  toward  it.  As  grid  points  vertically  migrate  toward  their 
target,  they  are  subjected  to  repelling  forces  from  grid  points  above  or  below.  This  device  keeps 
grid  points  in  a  column  from  fusing.  A  “pure”  isopycnic  model  may  be  viewed  as  a  hybrid  model 
in  which  the  repelling  force  is  set  to  zero. 

While  hybrid  coordinates  -  hybrid  in  the  narrow  BB81  sense  -  lost  some  of  their  appeal  in 
ocean  modeling  with  the  arrival  of  robust  methods  for  handling  fused  grid  points,  i.e.,  massless 
coordinate  layers,  the  idea  gained  acceptance  in  meteorology.  A  hybrid-isentropic  weather  pre¬ 
diction  model  based  on  the  BB81  concept  (Bleck  and  Benjamin,  1993)  has  been  in  use  by  the  US 
Weather  Service  since  the  mid-1990s  under  the  acronym  Rapid  Update  Cycle  (RUC). 

From  today’s  perspective,  adoption  of  hybrid  coordinates  in  isopycnic  ocean  modeling  is  ne¬ 
cessitated  not  so  much  by  the  outcropping  problem  but  by  the  degeneracy  of  isopycnic  coordinate 
representation  in  unstratified  or  convectively  unstable  water  columns.  In  this  paper  we  discuss 
some  early  results  obtained  with  a  hybrid  model  which  is  an  outgrowth  of  MICOM,  the  Miami 
Isopycnic  Coordinate  Ocean  Model  (Bleck  et  al.,  1992).  Emphasis  in  this  paper  is  on  a  fairly 
thorough  documentation  of  model  numerics;  hence,  there  is  no  room  here  for  an  exhaustive 
discussion  of  model  results.  Given  that  the  present  model  is  not  only  designed  for  specialized 
process  studies  but  also  for  comprehensive  tasks  such  as  climate  prediction,  emphasis  in  Section  5 
will  be  on  aspects  of  the  modeled  thermohaline-driven  circulation. 


2.  Model  overview 

In  this  section,  the  salient  features  of  the  new  hybrid  model  will  be  described.  The  model  will  be 
referred  to  as  HYCOM,  the  hybrid  coordinate  version  of  MICOM. 

2.1.  Governing  equations 

HYCOM,  like  MICOM,  is  a  primitive-equation  model  containing  5  prognostic  equations  -  two 
for  the  horizontal  velocity  components,  a  mass  continuity  or  layer  thickness  tendency  equa¬ 
tion,  and  two  conservation  equations  for  a  pair  of  thermodynamic  variables,  such  as  salt  and 
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temperature  or  salt  and  density.  The  model  equations,  written  in  (x,y,s)  coordinates,  where  5  is  an 
unspecified  vertical  coordinate,  are 
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where  v  =  ( u ,  r)  is  the  horizontal  velocity  vector,/?  is  pressure,  9  represents  any  one  of  the  model’s 
thermodynamic  variables,  a  =  p  '  is  the  potential  specific  volume,  £  =  ov/cx,  -  cu/cys  is  the 
relative  vorticity,  M  =  gz  +  pa  is  the  Montgomery  potential,  gz  =  cj>  is  the  geopotential,  /  is  the 
Coriolis  parameter,  k  is  the  vertical  unit  vector,  v  is  a  variable  eddy  viscosity/diffusivity  coefficient, 
and  t  is  the  wind-  and/or  bottom-drag  induced  shear  stress  vector.  represents  the  sum  of 
diabatic  source  terms,  including  diapycnal  mixing,  acting  on  6.  Subscripts  indicate  which  variable 
is  held  constant  during  partial  differentiation.  Distances  in  x,y  direction,  as  well  as  their  time 
derivatives  x  =  u  and  y  =  v,  are  measured  in  the  projection  onto  a  horizontal  plane.  This  con¬ 
vention  renders  the  coordinate  system  nonorthogonal  in  3-D  space  but  eliminates  metric  terms 
related  to  the  slope  of  the  s  surface  (Bleck,  1978a). 

Other  metric  terms,  created  when  vector  products  involving  (V-)  or  (Vx)  are  evaluated  on  a 
nonCartesian  grid  (e.g.,  in  spherical  coordinates),  are  absorbed  into  the  primary  terms  by 
evaluating  vorticity  and  horizontal  flux  divergences  in  ( 1)— (3)  as  line  integrals  around  individual 
grid  boxes.  Note  that  applying  V  to  a  scalar,  such  as  v2/2  in  (1),  does  not  give  rise  to  metric 
terms. 

After  vertical  integration  over  a  coordinate  layer  bounded  by  two  surfaces  stopAbot  (the  ocean 
surface,  the  sea  floor,  and  all  interior  layer  interfaces  are  s  surfaces  in  this  context),  the  continuity 
equation  (2)  becomes  a  prognostic  equation  for  the  layer  weight  per  unit  area,  A p  =  /?bot  —  Pt op: 
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The  expression  (sdp/ds)  represents  the  vertical  mass  flux  across  an  s  surface,  taken  to  be  positive  if 
in  T/?  (downward)  direction. 

Multiplication  of  (1)  by  dp/ds  and  integration  over  the  interval  (stoPAbot),  followed  by  division 
by  Ap/ As,  changes  the  shear  stress  term  in  that  equation  into 


g  , 
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while  the  lateral  momentum  mixing  term  integrates  to 
(A pVlVs  •  (vApV,v). 


(5) 


All  other  terms  in  (1)  retain  their  formal  appearance. 
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Wind-induced  and  bottom  stresses  are  prorated  among  coordinate  layers  in  accordance  with 
the  assumption  that  they  vary  linearly  over  a  prescribed  finite  depth  range  typically  chosen  to  be 
of  order  10  m. 

The  layer-integrated  form  of  (3)  is 

^6 A p  +  Vs  •  (vOAp)  +  (^s ^ O^J  -  =  Vj  ■  (vA pVs0)  +  (6) 


The  above  prognostic  equations  are  complemented  by  several  diagnostic  equations,  including 
•  the  hydrostatic  equation 


0M 

an  equation  of  state  linking  potential  temperature  T,  salinity  S,  and  pressure  p  to  or1 
and 


•  an  equation  prescribing  the  vertical  mass  flux  (sdp/ds)  through  an  s  surface. 

The  last-mentioned  equation  controls  both  spacing  and  movement  of  layer  interfaces  and  thus 
comprises  the  essence  of  hybrid  coordinate  modeling;  we  refer  to  the  algorithm  built  around  this 
equation  as  the  “grid  generator”,  to  be  described  in  detail  in  Section  3. 

The  horizontal  pressure  gradient  term  (Vv/V/  —  pS7su.  =  aX7:p)  in  (1)  must  be  formulated  so  as  to 
properly  transmit  interface  pressure  torques  up  and  down  the  water  column.  (Not  adhering  to  this 
principle  would,  for  example,  cause  the  model  to  violate  the  Sverdrup  relation  between  wind  stress 
curl  and  meridional  barotropic  flow.)  The  scheme  developed  for  this  purpose  is  discussed  in 
Appendix  A. 

The  use  of  potential  specific  volume  in  (7)  and  in  the  definition  of  M  can  be  justified  by  ar¬ 
guments  similar  to  those  given  in  Sun  et  al.  (1999,  Section  4b). 


2.2.  Transport  and  mixing  processes 

Owing  to  the  fact  that  coordinate  layers  in  HYCOM  are  predominantly  isopycnic,  the  model  in 
its  present  implementation  shares  with  MICOM  several  algorithms  for  handling  transport  and 
mixing  in  3-D  model  space. 

The  prognostic  equations  are  time-integrated  using  the  split-explicit  treatment  of  barotropic 
and  baroclinic  modes  developed  for  MICOM  (Bleck  and  Smith,  1990,  hereafter  BS90).  The  split- 
explicit  approach,  while  fraught  with  numerical  stability  problems  (Higdon  and  Bennett,  1997), 
has  proven  to  be  advantageous  for  executing  ocean  models  on  massively  parallel  computers  be¬ 
cause  it  does  not  require  solution  of  an  elliptic  equation. 


2.2.1.  Horizontal  transport  and  mixing 

Following  MICOM  tradition,  horizontal  mass  fluxes  are  computed  using  the  Flux  Corrected 
Transport  scheme  (Zalesak,  1979)  while  horizontal  tracer  transport,  which  in  a  model  featuring 
variable  vertical  mesh  spacing  must  be  done  in  flux  form,  is  handled  by  a  variant  of  the  MPDATA 
scheme  described  by  Drange  and  Bleck  (1997).  The  Coriolis  and  horizontal  momentum  advection 
terms  combined  [terms  2  and  3  in  (1)]  are  evaluated  in  potential-enstrophy  conserving  form.  To 
achieve  this,  the  finite-difference  analog  of  the  term  (£  +  /) k  x  v  must  be  written  as  an  analog  of 
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(((  +  f)/Ap)k  x  (vAp).  The  method  developed  for  handling  a  vanishing  denominator  (A p  =  0)  is 
described  in  Appendix  B  of  BS90. 

In  contrast  to  recent  MICOM  versions  where,  for  the  sake  of  computational  efficiency,  the 
constraint  ppot  =  const,  is  used  to  reduce  the  number  of  thermodynamic  prognostic  variables  from 
two  to  one  in  interior  layers,  HYCOM  carries  two  prognostic  thermodynamic  variables  every¬ 
where.  (The  prognostic  variable  used  in  MICOM  is  salinity.  Temperature  is  computed  by  in¬ 
verting  the  equation  of  state  which  for  this  purpose  is  approximated  by  a  third  degree  polynomial 
in  T  (Brydon  et  al.,  1999).) 

Intralayer  subgridscale  turbulent  mass  redistribution,  often  referred  to  as  bolus  transport,  is 
implemented  in  both  models  via  layer  interface  smoothing.  This  scheme  is  conceptually  identical 
to  the  Gent  and  McWilliams  (1990)  mixing  scheme,  except  that  it  is  used  here  in  its  underlying 
isopycnal  form.  Since  interface  smoothing  has  no  effect  in  regions  where  coordinate  surfaces  are 
geometrically  constrained,  the  bolus  transport  parameterization  in  HYCOM  is  only  active  in  the 
isopycnic  part  of  model  space.  This  restriction  is  unlikely  to  have  adverse  effects  on  model  dy¬ 
namics  as  long  as  nonisopycnic  coordinate  layers  are  essentially  confined  to  the  surface  mixed 
layer. 

The  interface  smoothing  term  in  MICOMs  continuity  equation  was  originally  formulated  as  a 
Laplace  operator,  meaning  that  interface  displacement  was  set  proportional  to  the  Laplacian  V2p 
of  interface  pressure  p.  Several  years  ago,  a  switch  was  made  to  biharmonic  smoothing  where 
interface  displacement  is  set  proportional  to  —V4p.  This  change  was  prompted  by  concerns  that 
Laplacian  interface  smoothing  may  have  a  systematic  shoaling  effect  on  layer  interfaces,  con¬ 
sidering  the  bowl  shape  of  isopycnals  in  the  vertical-meridional  plane. 

The  biharmonic  smoothing  concept  is  retained  in  HYCOM. 

2.2.2.  Surface  mixed  layer 

The  assured  presence  of  nicely  spaced  coordinate  layers  in  the  uppermost  part  of  the  water 
column  allows  formulation  of  turbulent  near-surface  mixing  in  terms  of,  for  example,  K  theory. 
For  reasons  of  computational  economy  and  model-to-model  continuity,  the  HYCOM  version 
discussed  here  retains  a  Kraus  and  Turner  (1967)  slab  mixed  layer. 

The  Kraus-Turner  (KT  hereafter)  closure  scheme  is  particularly  well-suited  to  models  that  treat 
the  mixed  layer  bottom  as  a  coordinate  surface.  This  is  the  case  in  MICOM  but  certainly  not  in 
the  present  model.  Attempts  to  suppress  the  diffusive  effect  of  the  spatial  mismatch  between  model 
interfaces  and  the  mixed  layer  bottom  -  mainly  by  splitting  the  coordinate  layer  bracketing  the 
mixed  layer  base  into  two  sublayers  -  were  altogether  unsuccessful.  The  KT  implementation 
chosen  for  the  present  suite  of  experiments  therefore  makes  no  effort  in  this  direction,  nor  does  it 
keep  a  record  of  the  depth  of  the  physical  mixed  layer  base  or  the  density  contrast  across  it  from 
one  time  step  to  the  next.  Rather,  it  treats  the  ad  hoc  stairstep  density  profile  defined  by  the 
discrete  model  layers  at  any  instant  as  “reality”  and  operates  on  it  accordingly.  (Note  that  so- 
called  layer  models  typically  treat  coordinate  layers  as  vertically  homogeneous;  hence,  properties 
such  as  temperature  and  salinity  are  piecewise  constant  in  the  vertical  but  exhibit  zero-order 
discontinuities  at  layer  interfaces.) 

The  KT  closure  scheme  operates  in  one  of  two  modes  depending  on  the  sign  of  the  surface 
buoyancy  flux.  A  buoyancy  flux  that  stabilizes  the  water  column  acts  as  a  drain  on  near-surface 
turbulence  kinetic  energy  (TKE).  In  this  case,  the  mixed  layer  depth  is  determined  by  the  ratio  of 
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the  TKE-producing  wind  stress  and  the  TKE-draining  buoyancy  flux,  that  is,  by  the  Monin  - 
Obukhov  length.  If  the  buoyancy  flux  destabilizes  the  water  column,  the  KT  scheme  stipulates 
that  the  TKE  generated  jointly  by  wind  stress  and  buoyancy  loss  is  converted  to  potential  energy 
by  entraining  into  the  mixed  layer  an  appropriate  amount  of  denser  water  from  below.  In  the  first 
case,  knowledge  of  the  mixed  layer  depth  at  the  beginning  of  a  given  model  time  step  is  not  re¬ 
quired,  but  in  the  second  case  the  outcome  of  the  entrainment  calculation  does  depend  on  this 
quantity. 

This  is  where  the  present  KT  implementation  sacrifices  accuracy  for  simplicity’s  sake.  Starting 
at  the  top  of  the  water  column,  the  program  looks  for  the  first  “stairstep”  in  the  density  profile 
and  defines  the  interface  where  this  happens  as  the  mixed  layer  base.  A  possible  concern  is  the  ill- 
posedness  introduced  by  the  fact  that  infinitesimal  changes  in  layer  density  may  lead  to  potentially 
large  changes  in  the  mixed  layer  base  estimate.  Fortunately,  this  ill-posedness  is  inconsequential  in 
the  context  of  the  KT  scheme  because  the  potential  energy  gained  by  entraining  water  of  density 
p  +  e  into  a  layer  of  density  p  approaches  zero  as  e  goes  to  zero.  This  is  to  say  that  the  outcome  of 
the  entrainment  calculation  is  not  overly  sensitive  to  misdiagnosing  the  initial  mixed  layer  base. 

Once  the  mixed  layer  depth  at  the  end  of  a  given  model  time  step  has  been  determined,  the 
water  column  is  homogenized  down  to  that  depth.  The  resulting  T /S  profile  is  then  projected  back 
onto  the  model  coordinate  layers  which  maintain  their  thickness  throughout  the  process.  This  last 
step,  which  is  not  found  in  MICOM,  generally  causes  water  mass  properties  to  be  exchanged 
across  the  mixed  layer  base.  It  makes  the  KT  mixed  layer  in  the  present  model  inherently  more 
diffusive  than  in  MICOM  and  calls  into  question  the  fidelity  of  the  KT  scheme  in  coarse-vertical 
resolution  models  unless  they  treat  -  as  MICOM  does  -  the  mixed  layer  base  as  a  layer  interface. 
In  an  attempt  to  counteract  this  diffusion,  the  TKE  generation  term  in  the  present  suite  of  ex¬ 
periments  has  been  artificially  reduced  by  50%. 


2.2.3.  Interior  diapycnal  mixing 

MICOM  uses  the  diapycnal  mixing  algorithm  of  McDougall  and  Dewar  (1998)  which  is  tai¬ 
lored  to  layer  models  treating  density  as  piecewise  constant  in  the  vertical.  Since  the  latter  is  also 
true  in  HYCOM,  the  Dewar-McDougall  scheme  is  carried  over  without  modification. 

In  the  HYCOM  version  discussed  here,  the  diapycnal  mixing  equations  are  solved  explicitly  in 
time.  Numerical  stability  constraints  in  this  case  do  not  allow  mixing  coefficients  of  the  size 
typically  found  in  turbulent  flow  regions.  Newer  versions  of  HYCOM  are  available  which,  by 
solving  the  relevant  equations  implicitly  (Hallberg,  2000),  maintain  numerical  stability  and  thus 
enable  the  model  to  use  K  theory  to  mix  fluid  throughout  the  water  column  regardless  of  whether 
the  flow  is  turbulent  or  laminar.  Both  mixed  layer  closure  and  patchy  turbulence  in  sill  overflows 
stand  to  benefit  from  this  treatment. 


2.2.4.  Deep  convection 

MICOM  maintains  stable  stratification  in  each  grid  column  by  absorbing  into  the  mixed  layer 
any  mass-containing  interior  layer  whose  density  falls  below  the  mixed  layer  density.  This  type  of 
mixed-layer  entrainment,  akin  to  convective  adjustment,  is  independent  of  entrainment  resulting 
from  TKE  closure.  It  plays  a  vital  role  in  creating  ultra-deep  “mixed  layers”  which  the  model 
needs  to  form  and  detrain  abyssal  water  masses. 
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HYCOM,  likewise,  resorts  to  convective  adjustment  to  maintain  stable  stratification.  However, 
since  density  in  this  model  is  not  guaranteed  to  increase  with  depth  in  layers  below  the  surface 
layer,  convective  adjustment  must  be  allowed  to  take  place  between  any  two  adjoining  model 
layers.  As  in  fixed-grid  models,  interface  depths  remain  unchanged  in  HYCOM  while  homoge¬ 
nization  of  the  pair  of  statically  unstable  layers  takes  place. 

2.3.  Sea  ice  model 

A  simple  sea  ice  model  is  used  is  to  handle  the  energetic  consequences  of  water  phase  changes  in 
polar  regions.  The  model,  which  has  much  in  common  with  one  developed  by  Semtner  (1976, 
Appendix)  focuses  on  three  (and  only  three)  aspects,  namely  the  prevention  of  water  temperatures 
below  the  freezing  point,  the  effect  of  freezing  and  melting  on  mixed-layer  salinity,  and  the  impact 
of  the  ice  surface  on  ocean-atmosphere  energy  fluxes. 

If  outgoing  surface  heat  fluxes  cause  the  mixed  layer  temperature  at  a  given  location  to  drop 
below  the  freezing  point  of  - 1 .8  °C,  the  ice  model  extends  a  repayable  “energy  loan”  to  the  mixed 
layer  at  this  location  to  maintain  a  — 1.8  °C  temperature.  The  loan  amount  obviously  is  pro¬ 
portional  to  the  amount  of  liquid  water  converted  to  ice.  If  the  surface  heat  flux  is  directed  so  as  to 
supply  heat  to  the  ocean,  any  outstanding  energy  loan  balance  at  a  given  location  must  be  repaid 
(i.e.,  the  existing  ice  must  be  melted)  before  SST  is  allowed  to  rise  above  —1.8  °C. 

Salinity  fluxes  though  the  ice-water  interface  that  result  from  melting  and  freezing  of  ice  are 
computed  based  on  the  assumption  that  ice  contains  25  g/kg  less  salt  than  sea  water. 

Sea  ice  profoundly  affects  surface  energy  fluxes,  both  by  virtue  of  its  high  albedo  compared 
to  water  and  because  an  ice  surface  can  be  much  colder  than  open  water.  The  approach  taken 
to  compute  heat  fluxes  across  the  upper  and  lower  surfaces  of  the  ice  sheet  is  described  in 
Appendix  B. 


3.  The  grid  generator 

Here  we  describe  an  experimental  scheme  that  allows  coordinate  layers  to  maintain  finite 
thickness,  at  the  price  of  becoming  nonisopycnic,  wherever  lack  of  stratification  would  cause  is- 
opycnic  layers  to  collapse  to  zero  thickness.  1  The  resulting  hybrid  coordinate  system  puts  iso- 
pycnic  grid  points  which  in  a  “pure”  isopycnic  model  are  deactivated  (i.e.,  fused  to  others  at  the 
sea  surface  or  the  bottom  of  the  slab  mixed  layer)  to  use  in  improving  vertical  model  resolution. 

An  important  point  to  note  is  that  layer  inflation  only  takes  place  in  the  upper  ocean;  no  efforts 
are  made  to  inflate  massless  coordinate  layers  fused  to  the  sea  floor.  Doing  so  would  create  steeply 
inclined  nonisopycnic  coordinate  layers  in  places  like  the  shelf  break  which  would  introduce  a 
class  of  numerical  problems  familiar  from  scaled-depth  (“cr”)  coordinate  2  models  (Janjic,  1977). 


1  Both  MICOM  and  HYCOM  sidestep  the  Boussinesq  approximation  by  using  pressure  p  in  lieu  of  height  z  as  the 
“geometric”  vertical  coordinate.  Accordingly,  terms  like  Cartesian,  thickness,  depth,  etc.,  refer  to  p  rather  than  z. 

2  The  letter  a  has  three  different  meanings  in  this  paper:  sea  water  density  anomaly,  pressure  divided  by  bottom 
pressure,  and  Stefan-Boltzmann  constant. 
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MICOMs  technique  for  handling  the  intersection  of  isopycnals  with  topography,  documented  in 
detail  in  Appendix  A  of  BS90,  has  withstood  the  test  of  time  (e.g.,  Smith,  1992)  and  has  been 
retained  in  HYCOM. 

Aside  from  maintaining  minimum  layer  thickness  in  the  upper  part  of  the  water  column,  the 
grid  generator  is  also  in  charge  of  restoring  isopycnic  conditions  in  coordinate  layers  wherever 
possible.  The  biggest  challenge  in  this  regard  is  the  annual  buildup  and  obliteration  of  the  sea¬ 
sonal  thermocline.  Warming  and  cooling  of  the  ocean  from  above  creates  coordinate  layers  in 
HYCOM  that  are  isopycnal  in  character  during  summer  but  change  into  geometrically  defined 
layers  in  winter.  Considerable  freedom  exists  in  placing  the  latter  in  the  water  column.  Since 
seasonal  heating/cooling  creates/destroys  isopycnals  right  at  the  surface,  a  strategy  allowing 
geometrically  defined  layers  to  stay  near  the  surface  during  winter  minimizes  seasonal  grid  mi¬ 
gration  and  concomitant  dispersion  of  water  mass  properties.  On  the  other  hand,  allowing  them 
to  descend  to  greater  depths  in  winter  improves  vertical  grid  spacing  in  deep  mixed  layers  and 
prevents  “data  voids”  that  may  cause  the  model  to  misrepresent  springtime  mode  water  forma¬ 
tion.  The  scheme  decribed  here  sacrifices  uniform  vertical  resolution  in  the  mixed  layer  for  the 
sake  of  minimizing  nonmaterial  interface  migration  during  spring  and  fall. 

The  scheme,  outlined  in  greater  detail  in  Appendix  C,  is  designed  to  be  highly  local  for  com¬ 
putational  simplicity’s  sake.  Not  only  does  it  operate  in  one  dimension  -  the  vertical  but  it 
determines  the  need  for  moving  gridpoints  (“regridding”)  by  considering  the  properties  of  only 
three  adjoining  coordinate  layers.  It  begins  by  checking  the  density  in  a  particular  layer  against  its 
target  value.  If  a  discrepancy  is  found,  only  one  interface  is  moved,  resulting  in  mass  exchange 
with  either  the  layer  above  or  below. 

Conceptually  the  simplest  way  of  reducing/increasing  water  density  in  a  given  layer  is  by  di¬ 
luting  the  layer  with  water  from  above/below.  This  approach  is  akin  to  upstream  or  donor  cell 
advection  of  mass  in  the  vertical,  and  as  such  it  has  rather  benign  (though  diffusive)  numerical 
properties.  It  also  has  the  potential  disadvantage  of  always  causing  layers  whose  density  is  off 
target  to  grow  at  the  expense  of  neighboring  layers. 

An  alternative  approach,  this  one  retaining  the  flavor  of  centered  advection  schemes,  is 
(i)  to  divide  the  layer  whose  density  we  wish  to  modify  into  two  sublayers  of  different  densities, 
and  (ii)  to  expel  one  sublayer.  For  this  scheme  to  yield  the  desired  result,  one  sublayer  density 
must  be  chosen  to  correspond  to  the  desired  target  density,  while  the  density  of  the  other  sublayer 
should  match  the  density  (actual  or  target)  of  the  layer  above  or  below,  as  the  case  may  be. 

While  it  may  sound  physically  implausible  that  the  density  in  a  layer  can  be  manipulated  by 
extruding  water  into  an  adjacent  layer,  this  operation  is  legitimate  as  long  as  it  does  not  create  new 
extrema  in  temperature  or  salinity.  For  this  reason,  the  “unmixing”  process  just  described  must 
remain  confined  to  a  “bounding  box”  in  T/S  space  whose  four  sides  are  defined  in  terms  of  the 
T/S  values  in  the  given  layer  and  its  two  neighbors.  See  Appendix  C  for  numerical  details. 

The  unmixing  (“UX”)  approach  just  described  differs  from  the  diluting  (“DL”)  approach 
described  earlier  in  that  it  causes  layers  whose  density  is  off  target  to  shrink.  Having  two  schemes 
available  allows  us  to  design  composite  regridding  schemes.  Given  that  in  many  situations  the  T/S 
bounding  box  constraint  will  prevent  formation  of  sublayers  having  the  appropriate  densities,  it  is 
necessary  to  follow  the  UX  step,  if  employed  at  all,  by  a  DL  step. 

The  effect  of  the  vertical  transport  terms  on  the  prognostic  variables  in  (4)  and  (6)  is  a  by¬ 
product  of  the  regridding  process  just  described.  Vertical  momentum  transport  is  evaluated 
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separately,  using  an  upstream  (donor  cell)  scheme  and  the  mass  flux  (sdp/ds)  implied  by  re- 
gridding. 

The  hybrid  grid  treatment  described  here  may  be  viewed  as  a  special  1-D  case  of  the  Arbitrary 
Lagrangian-Eulerian  (ALE)  scheme  introduced  by  Hirt  et  al.  (1974).  It  differs  from  common 
applications  of  the  ALE  scheme  in  one  important  aspect,  namely,  the  emphasis  placed  here  on 
regridding  for  the  sake  of  restoring  isopycnal  conditions  in  coordinate  layers,  in  addition  to  re- 
gridding  for  the  sake  of  maintaining  minimum  grid  point  spacing. 


4.  Experimental  setup 

HYCOMs  performance  will  be  demonstrated  here  by  a  series  of  100-year  simulations  carried 
out  on  a  near-global  domain  extending  from  69°S  to  65°N  with  a  horizontal  mesh  size  of 
1.4°  x  1.4°  cos  {latitude).  Water  depth  is  based  on  ET0P05  data,  averaged  over  individual  grid 
squares  but  otherwise  unsmoothed. 

The  vertical  structure  of  the  ocean  is  represented  by  14  layers  whose  target  densities  are  given  in 
Table  1.  Results  will  be  compared  with  those  from  MICOM  configured  for  the  same  domain  and 
same  14  layers.  (MICOMs  top  layer  is  a  variable-density  mixed  layer,  but  coordinate  values  for 
the  isopycnal  layers  match  values  2-14  from  Table  1.)  Thermobaric  effects  on  water  density  (Sun 
et  al.,  1999)  are  incorporated.  The  diapycnal  mixing  coefficient  is  set  to  (2  x  10  7  m2  s  2)/N  where 
N  is  the  buoyancy  frequency.  Isopycnal  diffusivity  and  viscosity  values,  including  the  one  used  for 
thickness  diffusion  (interface  smoothing),  are  formulated  as  udAx  where  Ax  is  the  local  horizontal 
mesh  size  and  ud  is  of  order  0.01  m  s-1.  In  regions  of  large  shear,  isopycnal  viscosity  is  set  pro¬ 
portional  to  the  product  of  mesh-size  squared  and  total  deformation  (Smagorinsky,  1963),  the 
proportionality  factor  being  0.2. 

Surface  forcing  is  based  on  monthly  climatologies  of  the  following  fields: 

•  surface  air  temperature,  relative  humidity,  and  rms  wind  speed  from  GOADS  (Woodruff  et  al., 
1987); 

•  wind  stress  components  from  the  European  Centre  for  Medium-Range  Weather  Forecasting; 

•  net  radiation  from  Oberhuber  (1988)  atlas; 

•  precipitation  based  on  satellite  microwave  measurements  (Spencer,  1993). 

Instantaneous  forcing  values  are  obtained  from  the  monthly  fields  using  quasi-Hermite  inter¬ 
polation  in  time. 

Annually  averaged  runoff  from  14  major  rivers  totaling  0.43  Sv  is  added  as  point  sources  to  the 
precipitation  field  at  the  appropriate  river  outflow  locations.  Another  0.27  Sv  of  freshwater  input 
representing  Antarctic  ice  melt  is  distributed  uniformly  along  the  southern  domain  boundary. 

Table  1 

Target  densities  (a2  =  ppot  —  1000  kg  m-3)  of  the  14  model  layers 
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Except  for  wind  stress,  the  omission  of  freshwater  input  along  the  northern  boundary,  and  the 
absence  of  imposed  diapycnal  fluxes  mimicking  the  effect  of  bottom  water  formation  in  polar 
basins  outside  the  computational  domain,  the  above  forcing  fields  agree  with  those  used  in  pre¬ 
vious  near-global  MICOM  simulations  (e.g.,  Bleck,  1998). 


5.  Results 

We  begin  by  showing  in  Fig.  1  an  example  of  HYCOMs  coordinate  layer  structure  in  relation 
to  the  density  field.  The  feature  to  note  is  the  approximate  alignment  of  coordinate  layers  and 
isopycnals  in  the  stratified  portion  of  the  ocean  and,  in  contrast  to  this,  their  tendency  to  un¬ 
couple  near  the  surface  where  coordinate  layers  flatten  out  while  isopycnals  rise  to  the  sea 
surface.  (Density  is  actually  piecewise  constant  in  the  vertical  within  a  coordinate  layer,  so  some 
artistic  license  has  been  taken  in  generating  this  plot.  The  apparent  density  inversions  in  the 
transition  zone  between  the  isopycnal  and  constant-depth  coordinate  domains,  as  well  as  in¬ 
stances  of  obvious  misalignment  of  the  two  sets  of  curves  in  the  deep  ocean,  are  caused  by  the 
vertical  spline  fitting  scheme  employed.  In  other  words,  in  locations  where  layer  interfaces  and 
density  isopleths  appear  to  be  only  approximately  aligned,  they  are  actually  perfectly  aligned  in 
the  model.) 

Fig.  1  illustrates  two  additional  points:  (i)  the  thickness  of  the  top  layer  is  held  constant,  except 
in  very  shallow  water  where  interfaces  are  compressed  in  a  coordinate  fashion,  and  (ii)  no  attempt 
is  presently  made  to  widen  the  vertical  spacing  of  coordinate  layers  in  unstratified  parts  of  the 
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Fig.  1.  Vertical  section  through  upper  1000  m  of  DL  solution  at  57.7°W.  Shaded  contours:  density  anomaly  (ay  units). 
Solid  lines:  coordinate  surfaces.  Small  numbers:  layer  indices.  Depth  marked  at  right,  latitude  along  bottom. 
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ocean,  such  as  the  Labrador  Sea  shown  on  the  left.  One  reason  for  adhering  to  minimal  spacing 
regardless  of  ambient  stratification  was  mentioned  in  Section  3.  Another  reason  is  discussed  in 
Appendix  C. 

Inspection  of  numerous  cross-sections  of  the  type  shown  in  Fig.  1  has  failed  to  reveal  any  signs 
of  layer  structure  deterioration  during  the  100-year  integration  period.  Specifically,  the  present 
grid  generator  shows  no  tendency  to  “coarsen”  the  steps  in  the  density  profile  by,  for  example, 
deflating  alternate  coordinate  layers.  Note  that  there  are  no  a-priori  safeguards  in  the  present  grid 
generator  that  would  prevent  degeneracies  and  hysteresis  effects  of  various  kinds  in  the  placement 
of  isopycnal  layers  in  the  water  column.  The  fact  that  the  model  is  capable  of  maintaining  a 
physically  consistent  layer  structure  over  long  periods,  though  difficult  to  document,  is  perhaps 
our  most  significant  finding. 

As  mentioned  earlier,  two  coordinate  maintenance  or  grid  generator  schemes,  named  DL  and 
UX,  are  being  compared.  The  DL  scheme  modifies  the  density  in  a  layer  by  diluting  it  with  water 
from  an  adjacent  layer,  while  the  UX  scheme  resorts  to  unmixing,  i.e.,  the  formation  of  two 
sublayers  of  different  densities  and  subsequent  ejection  of  one  sublayer,  to  modify  the  density  in  a 
given  layer. 

In  addition,  two  versions  of  HYCOM  are  compared,  one  solving  the  prognostic  equation  (3) 
for  potential  temperature  T  and  salinity  S,  and  one  solving  (3)  for  potential  density  ppot  and  S. 
The  two  model  versions  will  be  referred  to  as  TS  and  RHO,  respectively.  While  the  eddy  mixing 
term  on  the  right-hand  side  of  (3)  is  incorrect  if  ppot  replaces  T  as  dependent  variable,  treating  the 
pair  ppot,5  as  prognostic  variables  greatly  simplifies  the  task  of  the  grid  generator  whose  actions 
are  governed  by  the  vertical  distribution  of  ppot  rather  than  T  or  S.  As  pointed  out  in  Appendix  C, 
coordinate  maintenance  in  the  case  of  spatially  varying  T,S  is  inexact  as  the  amount  of  mass 
transferred  among  layers  is  determined  from  the  ppot  profile  while  the  properties  actually  trans¬ 
ferred  are  T  and  S. 

The  TS  and  RHO  versions  of  HYCOM  combined  with  the  DL  and  UX  grid  generator  options 
create  a  total  of  four  model  variants.  Space  restrictions  only  permit  presentation  of  a  few  key 
results  from  the  four  model  runs. 

A  relevant  question  to  ask  is  which  of  the  HYCOM  versions  tested  comes  closest  to  MI- 
COM  in  preserving  -  or  modifying,  as  the  case  may  be  -  the  oceanic  density  structure.  To  give 
an  impression  of  the  100-year  drift  in  the  MICOM  run  itself,  we  show  in  Fig.  2  the  contrast 
between  zonally  averaged  interface  depths  at  year  1  and  year  100.  The  model  drift  indicated 
in  the  figure  -  rising  interfaces  (cooling)  at  thermocline  depths  combined  with  interface  sink¬ 
ing  (warming)  in  the  abyssal  ocean  -  is  typical  of  century-scale  simulations  in  which  the 
diapycnal  mixing  coefficient  is  set  inversely  proportional  to  the  buoyancy  frequency.  Insuffi¬ 
cient  Antarctic  bottom  water  production  clearly  aggravates  the  density  drift  in  this  particular 
model  run. 

In  assessing  the  diffusive  properties  of  HYCOM  one  needs  to  keep  in  mind  that  coordinate 
maintenance  leads  to  nonzero  vertical  transport  terms  [the  terms  involving  s  in  (1)— (6)],  and  that 
one  of  the  grid  generators  evaluates  these  terms  using  the  notoriously  diffusive  upstream  or  donor 
cell  scheme.  Hence,  we  may  expect  the  RHO  version  of  HYCOM,  in  which  coordinate  mainte¬ 
nance  invoked  by  density  drift  due  to  cabbeling  does  not  occur,  to  be  less  diffusive  than  the  TS 
version,  at  least  in  the  deep  ocean.  Furthermore,  the  UX  grid  generator  scheme  can  be  expected  to 
be  less  diffusive  than  the  DL  scheme,  because  the  concept  of  “restoration  by  dilution”  itself  has  a 
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zonal  average  year  100.00  (jan.15)  MIC.002/N  <ref:  year  1.00  MIC.002/N> 


Fig.  2.  Vertical-meridional  section  showing  zonally  averaged  layer  interface  depths  in  MICOM  at  year  100.  Light/dark 
shading  indicates  downward/upward  interface  displacement  relative  to  interface  position  at  year  1.  Depth  marked  at 
right,  latitude  along  bottom. 


strong  donor  cell  flavor  whereas  the  concept  of  “restoration  by  unmixing”  does  not.  In  summary, 
the  TS/DL  combination  should  be  the  most  diffusive  one  while  the  RHO/UX  combination  should 
be  least  diffusive. 

The  above  expectations  turn  out  to  be  true.  To  illustrate  this,  we  show  in  Figs.  3  and  4  the 
zonally  averaged  vertical  density  structure  at  year  100  from  the  two  runs  showing  the  most  ex¬ 
treme  behavior,  TS/DL  and  RHO/UX,  referenced  to  the  MICOM  density  field  at  year  100.  To 
make  this  comparison  meaningful,  both  MICOM  and  HYCOM  results  have  been  converted  to 
potential  density  space  by  a  coordinate  transform  method  outlined  in  Appendix  D.  In  the  case  of 
MICOM,  this  transformation  boils  down  to  apportioning  the  nonisopycnic  uppermost  model 
layer  among  the  isopycnic  ones. 

Figs.  3  and  4  indicate  warming  in  most  parts  of  the  basin  compared  to  the  MICOM  run,  but  the 
warming  is  noticeably  stronger  in  the  TS/DL  run  than  the  RHO/UX  run.  In  both  HYCOM 
versions,  the  incremental  drift  appears  to  be  comparable  in  magnitude  to  the  drift  seen  in  MI¬ 
COM  itself.  We  interpret  this  as  an  indication  that  none  of  the  four  coordinate  maintenance 
schemes  tested  does  inordinate  damage  to  the  oceanic  density  structure. 

The  100-year  integration  period  is  long  enough  to  allow  for  a  significant  deterioration  of  the 
thermohaline-forced  meridional  overturning  circulation  (MOC)  in  an  ocean  model.  To  investigate 
this  aspect  of  model  performance,  we  plot  in  Fig.  5  stream  functions  of  the  zonally  integrated 
MOC  from  the  MICOM  run  in  the  three  major  ocean  basins  and  the  world  ocean  as  a  whole, 
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Fig.  3.  Vertical-meridional  section  showing  zonally  averaged  layer  interface  depths  in  TS/DL  run  at  year  100.  Light/ 
dark  shading  indicates  down  ward/up  ward  interface  displacement  relative  to  MICOM  reference  solution.  Depth  marked 
at  right,  latitude  along  bottom. 


averaged  over  model  years  100-110.  Fig.  6  shows  the  corresponding  results  from  the  HYCOM 
version  that  has  just  been  shown  to  deviate  most  strongly  from  MICOM,  namely,  the  TS/DL 
version. 

The  time-integrated  mass  fluxes,  on  which  the  overturning  stream  functions  are  based,  were 
transformed  from  hybrid  to  ppot  space,  using  once  again  the  scheme  presented  in  Appendix  D. 
The  apparent  compression  of  the  streamfunction  isopleths  along  the  upper  edge  in  the 
Pacific  and  Global  panels  could  have  been  avoided  by  widening  the  range  of  output  densities 
in  the  coordinate  transform.  Unfortunately,  since  the  transform  is  done  online,  i.e.,  during 
formation  of  the  time  integral,  repairing  this  defect  would  have  required  repeating  all  model 
runs. 

Several  features  in  Figs.  5  and  6  are  in  at  least  qualitative  agreement  with  observational  esti¬ 
mates  of  the  MOC  in  the  various  parts  of  the  world  ocean  (e.g.,  Schmitz,  1996a, b).  Significant 
amounts  of  deep  water  are  formed  both  in  the  North  Atlantic  and  the  Southern  Ocean  while  the 
Pacific  MOC  shows  relatively  shallow  cells  driven  by  equatorial  upwelling.  The  overturning  cell 
driven  by  sinking  in  the  North  Atlantic  extends  far  into  the  southern  hemisphere.  As  already 
mentioned  while  discussing  Fig.  2,  far  too  little  Antarctic  Bottom  Water  is  produced  in  the  present 
model  runs,  and  the  small  amount  formed  is  not  able  to  penetrate  the  northern  hemisphere  as 
observed. 

The  northern-hemispheric  meridional  heat  flux  is  seen  to  be  significantly  higher  in  TS/DL 
(Fig.  6)  than  MICOM  (Fig.  5).  In  the  Atlantic,  one  reason  for  this  difference  is  the  behavior 
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zonal  average  year  100.00  (jan.15)  HYC/RHO/UX  <ref:  year  100.00  MIC.002/N> 


Fig.  4.  As  in  Fig.  3,  but  for  RHO/UX  run. 


of  the  North  Atlantic  Drift  Current  in  the  two  models.  Fig.  7  indicates  that  the  NADC  in 
HYCOM  takes  a  more  northerly  path  compared  to  MICOM,  lowering  the  surface  density  south 
of  Greenland  by  as  much  as  0.55  kg  m  3.  The  corresponding  temperature  increase  (not  shown)  is 
as  high  as  3  °C,  leading  to  a  marked  additional  heat  loss  to  the  atmosphere.  No  explanation  is 
offered  at  this  time  for  the  discrepancy  regarding  NADC  positioning. 

One  of  MICOMs  traditional  strengths  is  its  ability  to  cleanly  delineate  geographic  regions 
where  the  conversion  of  surface  water  to  deep  water  takes  place.  It  is  of  interest  to  investigate 
whether  and  to  what  extent  HYCOM  retains  this  advantage.  The  following  results  are  based  on 
methods  developed  by  Sun  and  Bleck  (2000)  for  quantitatively  diagnosing  3-D  mass  fluxes  in 
(x,y,p  t)  space.  The  assertion  to  be  tested  is  that  the  RHO  version  of  HYCOM  does  a  better  job 
than  the  TS  version  in  delineating  the  downward  branch  of  the  MOC.  As  before,  we  use  the 
MICOM  solution  as  reference  and  base  the  discussion  on  model  output  transformed  to  ppot  space 
and  averaged  over  model  years  100-110. 

Fig.  8  shows  that  MICOM  injects  roughly  9  Sv  of  water  into  isopycnal  layer  12  (and  possibly 
beyond)  in  the  Labrador  Sea,  but  rather  insignificant  amounts  elsewhere  in  the  North  Atlantic.  It 
is  worth  pointing  out  that  that  these  diapycnal  fluxes  are  far  too  large  to  result  from  interior 
diapycnal  mixing;  rather,  they  are  caused  by  water  mass  transformation  in  the  surface  mixed  layer 
which  reaches  into  isopycnal  layer  12  during  at  least  part  of  the  year. 

The  corresponding  field  from  the  TS/DL  experiment  is  displayed  in  Fig.  9.  The  string  of  up-  and 
downwelling  cells  along  the  Gulf  Stream  extension  axis  is  the  item  of  concern.  As  anticipated, 
these  centers  are  largely  eliminated  if  ppot,S  replace  T.  S  as  prognostic  thermodynamic  variables 
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Fig.  5.  Top  4  panels:  meridional  overturning  streamfunction  (Sv)  in  the  3  major  basins  and  the  world  ocean  from 
MICOM.  Bottom:  meridional  heat  flux  (PW)  in  the  3  major  basins  (AT,  IN,  PA)  and  the  world  ocean  (GL). 
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Fig.  6.  As  in  Fig.  5,  but  for  TS/DL  run. 
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Fig.  7.  North  Atlantic  surface  density  difference  (rr2  units)  between  TS/DL  and  MICOM  run  at  year  100.  Latitude/ 
longitude  marked  at  left/bottom. 


(Fig.  10).  Note  the  substantial  downward  mass  flux  (11  Sv)  into  layer  12  in  the  Labrador  Sea  in 
the  RFIO/DL  simulation,  one  of  many  indicators  that  the  RHO  version  of  HYCOM  bears  closer 
resemblance  to  MICOM  than  does  the  TS  version. 

The  similarities  between  Figs.  8  and  10  are  no  accident  and  come  at  a  price.  Both  MICOM 
and  HYCOMs  RHO  version  achieve  their  relatively  clean  separation  of  3-D  mass  fluxes  into 
isopycnal  and  diapycnal  components  by  emphasizing  conservation  of  ppot  at  the  expense  of 
conserving  T.  Unfortunately,  this  strategy  eliminates  not  only  numerically-induced  cabbeling  - 
that  is,  cabbeling  caused  by  mutual  inconsistencies  in  T,S  advection  but  also  “physical” 
cabbeling  related  to  isopycnal  diffusion.  The  relative  importance  of  explicit  T,  S  conservation 
achieved  in  the  TS  version  of  HYCOM  and  the  clarity  of  the  model’s  rendering  of  the  ther- 
mohaline-forced  circulation  is  an  issue  whose  resolution  will  undoubtedly  depend  on  the  specific 
model  application. 

One  measure  of  success  of  an  oceanic  circulation  model  is  its  ability  to  portray  the  equatorial 
current  system.  The  performance  of  HYCOM  in  this  regard  is  illustrated  in  Fig.  11  where  me¬ 
ridional  cross-equatorial  sections  through  the  upper  300  m  of  the  water  column  from  the  TS/DL 
and  RHO/UX  runs  are  displayed  together  with  the  corresponding  section  from  MICOM.  As  in 
Fig.  1,  some  liberty  has  been  taken  in  drawing  contour  lines;  fields  that  are  piecewise  constant  in 
one  coordinate  direction  (the  vertical  in  this  case)  are  difficult  to  display  in  black-and-white  line 
graphics  mode. 
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diapyc.vel  (m/yr)  thru  botm.  of  layer  11,  year  100.00  —  110.00  MIC.002/N 


Fig.  8.  Diapycnal  flux  through  bottom  of  layer  1 1  in  MICOM  run.  Grey-shaded  contours:  diapycnal  vertical  velocity 
(m/yr),  positive  downward.  Solitary,  variable-size  numbers:  regionally  integrated  diapycnal  transport  (0.1  Sv).  Un¬ 
shaded  contour  lines:  mean  sea  surface  height  (cm).  Latitude/longitude  marked  at  left/bottom. 

The  feature  of  interest  is  the  equatorial  undercurrent,  i.e.,  the  closed  set  of  isotachs  centered  on 
the  equator  at  150  m  depth.  The  undercurrent  has  a  surprisingly  similar  appearance  in  the  three 
cross-sections,  but  fairly  large  differences  are  found  in  the  strength  of  the  equatorial  surface  flow. 
The  high  equatorial  surface  flow  speed  generated  by  MICOM,  associated  with  a  somewhat  ex¬ 
cessive  equatorial  SST  minimum  (note  the  strong  doming  of  layer  interfaces  just  below  the  surface 
in  the  MICOM  solution),  has  been  a  long-standing  concern.  Specific  reasons  for  this  behavior 
have  not  been  isolated.  Layer  interface  displacements  and  flow  speeds  in  the  HYCOM  solutions  in 
Fig.  11  create  the  impression  of  higher  implicit  diffusion  in  HYCOM  compared  to  MICOM. 
Qualitatively,  this  difference  is  in  accord  with  the  notion  that  HYCOM  behaves  like  a  fixed-grid 
model  at  shallow  depths. 


6.  Concluding  remarks 

An  attempt  is  made  to  combine,  within  a  single  model  framework,  certain  features  found  in 
fixed-grid  and  isopycnic  coordinate  ocean  circulation  models  that  are  usually  considered  mutually 
exclusive,  namely,  maintenance  of  vertical  grid  resolution  regardless  of  stratification  on  the  one 
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Fig.  9.  Diapycnal  flux  through  bottom  of  layer  11  in  TS/DL  run.  See  Fig.  8  for  details. 


hand,  and  alignment  of  coordinate  surfaces  with  isopycnals  on  the  other.  Emphasis  in  this  article 
is  on  technical  documentation,  but  results  from  a  few  tests  are  presented  that  address  the  model’s 
numerical  and  physical  resilience. 

The  most  unconventional,  perhaps  even  discomforting  aspect  of  the  present  model  is  the  lack  of 
a  concise  definition  of  the  vertical  coordinate.  HYCOM  grid  points  are  free  to  move  up  and  down 
in  the  water  column,  subject  only  to  the  “mandate”  to  maximize  isopycnal  coordinate  repre¬ 
sentation  in  the  model  domain  while  at  the  same  time  avoiding  formation  of  massless  layers  near 
the  sea  surface.  It  is  unlikely  that  a  closed  formula  exists  that  can  match  the  present  algorithm  in 
its  ability  to  optimize  the  effect  of  these  two  constraints  in  a  noninterfering  manner. 

The  “discomfort”  most  likely  will  be  with  us  for  some  time,  and  since  there  is  no  formal  as¬ 
surance  or  proof  of  uniform  and  predictable  model  behavior  in  situations  ranging  from  typical  to 
extreme  -  the  type  of  assurance  taken  for  granted,  for  example,  in  the  case  of  scaled-depth  (cr) 
coordinates  popular  in  coastal  and  atmospheric  models  the  best  course  of  action  for  the  time 
being  is  to  put  the  model  through  an  extensive  series  of  practical  tests.  The  results  presented  here 
should  be  viewed  as  a  first  round  of  such  tests. 

At  the  core  of  HYCOM  is  the  so-called  grid  generator,  an  algorithm  which  moves  grid  points 
vertically  through  the  fluid  if  it  senses  opportunities  for  improving  isopycnal  alignment  of  co¬ 
ordinate  surfaces,  or  sees  the  need  to  restore  their  minimum  separation.  While  the  particular  grid 
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Fig.  10.  Diapycnal  flux  through  bottom  of  layer  11  in  RHO/DL  run.  See  Fig.  8  for  details. 


generator  tested  here  appears  to  produce  acceptable  results  in  a  broad  sense,  it  should  not  be 
taken  as  the  final  word  with  regard  to  both  coordinate  surface  placement  and  handling  of  the 
vertical  advection  terms  resulting  from  nonmaterial  movement  of  grid  points.  Our  ultimate  goal  is 
to  provide  the  user  with  a  set  of  algorithmic  and  parametric  grid  generator  options,  and  to  en¬ 
courage  independent  experimentation. 

Substantially  more  work  is  also  needed  (i)  to  determine  the  model’s  sensitivity  to  the  mini¬ 
mum  layer  thickness  parameter  <5  used  in  the  grid  generator,  and  (ii)  to  develop  strategies  for 
finding  the  optimal  number  of  low-density  model  layers.  A  large  S  combined  with  too  many 
model  layers  assigned  to  low  target  densities  (i.e.,  to  the  warm  water  sphere)  increases  the 
likelihood  that  the  Cartesian  portion  of  the  model  domain  will  encroach  upon  the  stratified 
region  below  the  mixed  layer  where  isopycnal  coordinate  representation  clearly  is  desirable.  The 
opposite  extreme,  a  very  small  <5  combined  with  too  few  low-density  layers,  is  likely  to  degrade 
modeling  of  surface  mixed  layer  processes.  Investigation  of  the  various  tradeoff's  sketched  here 
has  barely  begun. 

The  present  article  totally  sidesteps  the  question  of  how  modern  mixed  layer  closure  schemes 
fare  in  a  hybrid  coordinate  environment.  The  guaranteed  presence  of  nicely  separated  coordinate 
surfaces  in  the  oceanic  surface  layer  invites  replacement  of  MICOMs  Kraus-Turner  bulk  mixed 
layer  scheme  by  one  based  on  second-order  closure  (K  theory)  or  even  higher  order  closure 
schemes.  Experience  with  K  theory  schemes  will  be  described  in  subsequent  articles. 
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merid.sec.1 25.1  6w  year  100.00  (jan.15)  MIC.002/N 


merid.sec.125.16w  year  100.00  (jan.15)  HYCOM/DLUT 


Fig.  1 1 .  Cross-equatorial  section  through  uppermost  300  m  of  MICOM  (top),  TS/DL  (middle),  and  RHO/UX  solution 
(bottom)  near  125°W.  Shaded  contours:  zonal  velocity  (contour  interval:  10  cm/s).  Depth  marked  at  right,  latitude 
along  bottom.  Layers  marked  3, 4,  5, . . .  correspond  to  model  layers  1, 2,  3, . . .  given  in  Table  1.  Layers  marked  1  and  2, 
assigned  densities  27.79  and  30.07,  respectively,  were  added  during  postprocessing  (see  Appendix  D)  to  assure  correct 
rendering  of  low  density  water. 
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Appendix  A.  Formulation  of  the  horizontal  pressure  gradient  force  (PGF)  in  generalized  coordinates 

Both  MICOM  and  HYCOM  solve  the  horizontal  momentum  equations  in  advective  form. 
Care  must  be  exercised  in  this  case  to  retain  certain  conservation  properties  inherent  in  the  flux 
form  of  the  equations. 

In  a  hydrostatic  fluid  (6 <f>/ds  =  —a dp/ds),  the  layer-mass  weighted  horizontal  PGF  can  be 
expressed  as 

^[a.Vsp  +  Vs(j)\  =  (A-!) 

The  form  on  the  right-hand  side  of  (A.l)  has  well-known  implications  for  vortex  spinup/spin- 
down.  Specifically,  given  that  the  curl  of  this  expression  is 

0 

—  (V,  x  pVs(p), 

we  can  state  that  interface  pressure  torques  governing  vortex  spinup/spindown  in  individual  s 
coordinate  layers  have  the  form  (V.s  x  pVs(j)). 

Orderly  transmission  of  pressure  torques  up  and  down  the  water  column  is  important  for 
maintaining  Sverdrup  balance.  The  task  before  us,  therefore,  is  to  find  a  finite-difference  ex¬ 
pression  for  the  PGF  term  \uS7sp  +  V.v r/j ]  in  the  horizontal  momentum  equation  that  can,  after 
multiplication  by  layer  thickness,  be  transformed  by  finite-difference  operations  into  a  finite-dif¬ 
ference  analog  of  the  right-hand  side  of  (A.l). 

We  start  by  writing  the  x  component  of  the  last  term  in  (A.l)  in  the  simplest  possible  -  and 
hence  noncontroversial  -  form  5s(j?8x<j)),  where  the  S  and  overbar  operators  denote  2-point 
centered  differentiation  and  interpolation  (averaging),  respectively.  Finite-difference  product  dif¬ 
ferentiation  rules  allow  this  expression  to  be  expanded  as  follows: 

Ss(px8x<p)  =  (Sspx)Sx(/>s  +pxs8x5s(j) 

=  ( Sspx)Sx4 -  pxs5x(a5sp) 

=  ( 8spx)3x(f)S  -  8Jpsa8sp)  +  aSsp'dxps. 
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A  finite-difference  equation  matching  (A.l)  is  now  obtained  by  rearranging  terms  and  adding  an 
analogous  expression  for  the  y  component: 

'i6sp3xp'  +  ( Sspx)Sx(j)S  =  5x(psa8sp)  +  8s(px8x(f)), 
tx8spdyps  +  ( 8spv)dy(j)S  =  8y(psoc8sp )  +  8s(p’8y(f)). 

The  last  two  equations  state  that,  in  order  to  preserve  the  conservation  properties  expressed  by 
(A.l),  the  finite-difference  PGF  must  be  evaluated  in  the  form 


aVsP  +  Vs(p 


^  =r  <XxPs  +  <V/>S  ^ 

osp 


y-8sp  , 


V  <5, 


SypS  +  8y(f) 


(A.2) 


The  salient  result  of  the  above  analysis  is  that  writing  the  undifferentiated  factor  a  in  the  PGF 
formula  as  simply  of  or  of  can  lead  to  spurious  momentum  and  vorticity  generation.  To  avoid  this 
pitfall,  a  must  appear  in  the  PGF  formula  in  layer  thickness-weighted  form. 

For  use  in  isopycnal  or  quasi-isopycnal  models,  it  is  convenient  to  express  the  PGF  in  terms  of 
the  Montgomery  potential  M  =  <f>  +  pa.  The  proper  finite-difference  analog  of  M  in  a  staggered 
vertical  grid  (p  and  4>  carried  on  layer  interfaces,  a  carried  within  a  layer)  is 

M  =  <j>  +  aps , 


which  can  easily  be  shown  to  be  consistent  with  the  two  common  forms  of  the  hydrostatic 
equation, 

5s(j)  =  —ocdsp,  8SM  =  p8sa. 


(These  are  finite-difference  analogs  of  d(j)/dp  =  —a  and  dM/da  =  p.  respectively.) 

We  begin  the  task  of  converting  (A.2)  into  a  form  involving  M  by  writing  the  x  component  of 
(A.2)  as 


C^£rdxpF  +  dx<j)S  =  SXM  +  8xps  -  5x(a.ps)  .  (A.3) 

8sp  l  8sp 

Making  use  of  the  relation 

W  -IT  =  '-(i,rMSjs), 

where  8’x  represents  the  difference  between  two  neighboring  grid  points,  i.e.,  8’x  =  A x8x,  the  term  in 
square  brackets  in  (A.3)  can  be  expanded  into 


t-t  ( rJ-8sp  -  <xx8sp  )  8xps  -  psx8xa  = 


8sp 


1 


4  8sp 

1 

4  8spx  L 


{8'x8sp)  (8'xa)8xps  -  psx8x a 


(5'x8,p)(5'jf)  - 


8X  a. 
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The  above  expression  involves  a  total  of  four  p  points,  located  one  grid  distance  Ax  apart  on  two 
consecutive  s  surfaces.  Substantial  simplification  of  this  expression  is  possible  by  labeling  the  four 
points  as 


P\ 


Pi 


f  Ax  As  \ 

K=T+Y’s_yJ’ 

(  Ax  As  \ 

p*=p(*  +  Y’s  +  Y  ' 


With  a  modest  amount  of  arithmetic,  it  can  now  be  shown  that  the  term  in  square  brackets  in 
(A.3)  reduces  to 

PiPi-PiPA  s 
(j>4  -Pi)  +  ( Pi  ~P\) 

This  term,  which  in  combination  with  the  term  8XM  gives  the  PGF  in  x  direction,  is  the  sought- 
after  finite-difference  analog  of  —pc> a/6x  in 

aFJsp  +  Vscj)  =  V^M  —  pS7sct.. 

The  finite-difference  expression  for  the  PGF  in  y  direction  is  analogous. 


Appendix  B.  Vertical  heat  flux  calculation  in  the  presence  of  ice 


In  the  simple  “energy  loan”  ice  model,  the  ice  surface  temperature  calculation  is  based  on  the 
assumption  that  the  system  is  energetically  in  a  steady  state,  i.e.,  the  heat  flux  through  the  ice 
matches  the  atmospheric  heat  flux.  To  illustrate  this  approach,  we  write  the  atmospheric  heat  flux 
as  Fair  =  a(  T  —  Ta),  and  the  heat  flux  through  the  ice  as  Fice  =  b(Tw  —  T),  where  Ta,  Tw  denote  air 
and  water  temperature,  respectively,  while  T  is  the  ice  surface  temperature  -  the  unknown  in  the 
problem.  a,b  are  proportionality  factors. 

Given  Ta,  Tw  and  a  first-guess  value  for  T.:  we  seek  to  modify  T  such  that  the  difference  between 
Fee  and  Fair  is  eliminated,  or  at  least  reduced.  The  equation  describing  this  situation  is 

a{T  +  d  T  -  T.d)  =  b(Tw  -T-dT), 


which  yields 


dT 


aTa  +  bTw 
a  +  b 


(B.l) 


In  order  to  make  this  formula  applicable  in  situations  where  b]ur  is  a  mixture  of  sensible,  latent, 
and  radiative  heat  fluxes,  we  add  aT  -  aT  and  bT  —  bT  to  the  numerator  of  (B.l): 

^  _  aTa  —  aT  +  aT  +  bTw  —  bT  +  bT  ^ 
a  +  b 


After  substituting  the  original  definitions  of  Fail,Fice,  this  expression  reduces  to  the  rather  general 
formula 
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dT  =  Fi  ce  Fail 
a  +  b 

The  new  temperature  T  +  dT  must  not  be  allowed  to  exceed  0  °C  until  the  ice  is  melted  com¬ 
pletely. 

Practical  application  of  (B.2)  requires  knowledge  of  the  coefficients  a,b  which  represent  the 
derivatives  dFair/dT  and  dFice/dr,  respectively.  Guidance  on  the  magnitude  of  a  can  be  obtained 
from  the  conventional  heat  flux  bulk  formula.  It  suggests  a  =  ctpcvU  where  ct  is  a  nondimensional 
transfer  coefficient  (similar  to  the  drag  coefficient),  p  is  the  air  density,  cp  is  the  specific  heat  of  air 
at  constant  pressure,  and  U  is  the  wind  speed.  The  formula  for  radiative  energy  loss,  cT4,  suggests 
that  the  above  estimate  of  a  should  be  incremented  by  an  amount  of  order  4oT3.  A  reasonable 
choice  for  b  is  the  ratio  of  ice  thermal  conductivity  and  ice  thickness,  kice/Hici,. 

In  the  interest  of  computational  efficiency  in  coupled  climate  models,  we  wish  to  minimize 
information  exchange  with  the  atmosphere.  Toward  this  end,  we  make  a  independent  of  atmo¬ 
spheric  state  variables.  To  avoid  oscillatory  behavior  in  (B.2),  a  must  be  chosen  somewhat  larger 
than  “typical”  values  of  ctpcpU  +  4oT3;  in  other  words,  we  adopt  a  strategy  of  prudent  under¬ 
relaxation  of  T. 

Returning  briefly  to  the  energy  loan  concept,  we  finally  need  a  statement  relating  the  rate  of 
energy  borrowing  or  repaying  to  the  composite  atmospheric  heat  flux 


F  =  cTjce  +  (1  —  c)^opw>  (B.3) 

where  F0 pw  is  the  energy  flux  over  open  water  and  c  is  the  fractional  ice  coverage. 

We  assume  that  F  as  defined  in  (B.3)  is  the  energy  flux  felt  by  the  ocean  irrespective  of  the 
presence  of  ice.  In  other  words,  we  assume  that  the  energy  flux  between  atmosphere  and  ice,  Fice, 
equals  the  energy  flux  between  ice  and  ocean.  This  assumption  is  compatible  with  the  steady  state 
(zero  heat  flux  divergence)  assumption  made  in  deriving  (B.2). 


Appendix  C.  A  prototype  grid  generator 

Consider  three  consecutive  isopycnic  layers,  labeled  0,  1,  2,  in  a  stratified  water  column.  Let 
ak  >  a.k+ 1,  where  a  is  the  specific  volume.  Label  the  interfaces  such  that  interface  k  denotes  the 
upper  interface  of  layer  k. 

Suppose  oq  differs  from  its  layer  reference  value  oq.  We  seek  ways  to  re-discretize  the  water 
column  in  a  manner  that  preserves  the  overall  height  of  the  column,  represented  by  the  integral 
/  ad p,  while  changing  to  oq.  Conservation  of  the  integral,  which  is  important  for  preserving  the 
geostrophic  balance  throughout  the  column,  requires  that  alterations  of  a,  must  be  accompanied 
by  vertical  displacement  (in  p  space)  of  one  or  more  interfaces.  The  simplest  solution  is  to  move 
only  one  interface  and  to  select  this  interface  in  the  spirit  of  the  “donor  cell”  transport  scheme.  By 
this  we  mean  that  the  upper  interface  is  moved  upward,  i.e.,  layer  1  is  diluted  with  lighter  water 
from  layer  0,  if  layer  1  is  too  dense  (oq  <  oq).  Analogously,  the  lower  interface  is  moved  down¬ 
ward,  i.e.,  layer  1  is  diluted  with  denser  water  from  layer  2,  if  layer  1  is  too  light  (oq  >  oq).  This 
method  of  restoration  does  not  work  in  the  special  case  where  the  coordinate  layer  touching  the 
sea  floor  is  too  light,  and  like  all  upstream  advection  schemes  it  is  inherently  diffusive.  Methods 
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for  conserving  the  integral  J'  a  dp  by  extruding  water  to  a  neighboring  layer,  rather  than  entraining 
water,  will  be  discussed  later. 


Case  1  (oq  <  iq).  This  is  the  case  where  the  upper  interface  is  moved,  i.e.,  mass  is  exchanged 
between  layers  0  and  1 .  Conservation  of  f  a  dp  requires 

ixo(pi  -po)  +  <xi(p2  ~p  1)  =  ao(Pi  - Po )  +  ai(/>2  ~P\),  (C.l) 

where  px  is  the  pressure  of  the  upper  interface  after  re-discretization.  Solving  (C.l)  for  px  yields  the 
expression 


Pi  = 


Pi(ocq  -  «i)  +j?2(«l  -  «l) 
(Xq  —  QL\ 


(C.2) 


Note  that  the  weight  assigned  to  p2  is  negative,  indicating  that  for  large  discrepancies  between  oq 
and  oq  (C.2)  will  not  necessarily  yield  a  solution  px  >  p0.  To  assure  a  minimum  interface  spacing, 
we  therefore  replace  px  by 

px  =  max(p!,po  +  A0),  (C.3) 


where  Ao  is  a  minimum  pressure  difference.  Moving  the  interface  to  px  instead  of  px  means,  of 
course,  that  layer  1  will  end  up  with  a  specific  volume  5.\  different  from  the  reference  value  oq. 


Following  BB81,  A0  is  chosen  to  be  a  continuously  differentiable  function  which  for  large 
positive  arguments  A p  =px  —  p0  returns  the  argument  A p  (meaning  that  /),  =  px),  while  for  large 
negative  arguments  it  returns  a  small  constant  5  representing  the  lowest  permissible  layer  thick¬ 
ness.  Our  traditional  choice  for  this  so-called  “cushion”  function  (Bleck  and  Benjamin,  1993)  is 


A p 


A0  =  < 


1  (  Ap 

1  +  3(2«+1 


if 


if 


2<f  <4- 


(C.4) 


which  is  obtained  by  fitting  a  piece  of  a  parabola  in  (Ap,  A0)  space  between  the  two  straight  lines 
A0  =  Ap  and  A0  =  5.  Another  example  from  this  class  of  spliced-together  cushion  functions,  this 
one  acting  over  a  somewhat  wider  range  of  A p/8,  is 


Ap 


A0  =  \ 


4  f  Ap 

1  +  5  (  45  + 1 


A p  , 

lf  T>6’ 

if  -  4  <  ^  <  6, 
o 


(C.5) 


if -4  > 


Ap 


Note  that  the  above  cushion  functions,  whose  purpose  it  is  to  smooth  the  transition  between 
isopycnally  and  geometrically  constrained  segments  of  a  coordinate  surface,  make  no  reference  to 
neighboring  grid  points.  This  simplicity  lends  a  degree  of  robustness  to  our  hybrid  scheme  not 
necessarily  found  in  more  complicated  formulas. 

Two  cases  must  be  distinguished  in  computing  6q  from  (C.l): 


82 


R.  Bleck  /  Ocean  Modelling  37  (2002)  55-88 


Case  la  (px  <  p{).  In  this  case,  the  interface  ends  up  at  a  lower  pressure,  allowing  lighter  water  to 
raise  oq  to  a  value  a i  closer  -  if  not  equal  -  to  the  reference  value  ai.  The  value  a \  is  obtained  by 
substituting  tildes  for  carats  in  (C.l): 


Oil 


U\{pi  - Pi )  +  «o(pi  -Pi) 
P2  ~  Pi 


(C.6) 


Case  lb  (px  >  p{).  In  this  case,  layer  1  “donates”  water  to  layer  0;  hence,  it  is  a0,  rather  than  a u 
that  is  being  modified.  Conservation  of  fa  dp  then  leads  to 

~  _  cciipi  -  Po)  +  a0{pi  -  p0) 

CCq  ^  /  J 

P\  -Po 

Before  solving  (C.7),  px  must  be  modified,  if  necessary,  to  satisfy  px^pi  -  A2,  where  A2  is  a  layer 
thickness  threshold  presently  set  to  {p2  —  p\)/2. 


Case  2  (oq  >  ai).  This  is  the  situation  in  which  we  decide  to  move  the  lower  interface.  The 
equation  replacing  (C.l)  in  this  case  is 


ai(i?2  ~ P\)  +  a-iipi  - Pi )  =  op{p2  ~Pi)  +  ar(p3  ~P2)- 
It  calls  for  moving  the  pressure  of  the  lower  interface  from  p2  to 
Pi{&i  ~  ai)  +P2{ui  -  oc2) 

Pi  ~ 

a.\  —  a2 


(C.8) 

(C.9) 


The  concern  in  this  case  is  that  p2  may  get  too  close  to,  or  even  exceed,  p2.  We  therefore  replace  p2 
by 

p2  =  min {p2,p2  -  A3)  (C.10) 

with  A 3  currently  set  to  (p2  —  p2)/ 2.  Replacing  p2  by  p2  in  (C.8)  changes  a.\  into 

-i  =  ^fe-K)+«,fe-w)^  (C11) 

Pi  ~P\ 


Case  2a  (oq  >  ax  where  layer  1  is  the  bottom  layer).  Since  the  lower  interface  cannot  move  down  in 
this  case,  we  use  (C.l)  and  compute  an  upper  interface  movement  from  (C.2).  While  it  may  sound 
physically  implausible  that  the  density  in  a  layer  can  be  manipulated  by  extruding  water  into  an 
adjacent  layer,  this  operation  is  legitimate  if  one  conceptually  allows  density  to  vary  vertically 
within  a  coordinate  layer. 

Specifically,  as  long  as  a0  >  oq(>  oh),  layer  1  may  be  viewed  as  being  composed  of  two  sub¬ 
layers,  one  of  specific  volume  a0  and  one  of  specific  volume  oq .  Having  the  same  density  as  layer  0, 
the  lighter  sublayer  can  be  merged  with  layer  0.  The  denser  sublayer,  which  has  the  desired  density 
ax,  stays  in  layer  1. 

To  keep  px  from  exceeding  the  bottom  pressure  p2,  we  replace  px  by 
Pi  =  min (pi,p2  -  A2) 
with  A2  currently  set  to  [p2  —  p\)/2. 


(C.12) 
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The  above  algorithm  for  redistributing  mass  and  density  among  model  layers  needs  to  be 
amended  by  a  prescription  for  vertically  redistributing  thermodynamic  variables  other  than  mass 
and  density.  This  is  a  particularly  vexing  issue  for  hybrid  models  carrying  T,  S  (as  opposed  to  a,  S ) 
as  prognostic  variables. 

One  possible  approach  is  to  translate  pk  —  pk  into  a  generalized  vertical  velocity  (i.e.,  vertical 
motion  relative  to  the  layer  interface)  which  is  then  used  to  advect  T  and  S.  Ideally,  this  vertical 
advection  process  should  (i)  conserve  T,  S  in  the  column,  and  (ii)  produce  T,  S  changes  in  layer  1 
that  satisfy  (C.l). 

A  scheme  that  exactly  satisfies  (i)  while  approximately  satisfying  (ii)  can  be  obtained  by  re¬ 
apportioning  T,  S  among  the  layers  in  a  manner  analogous  to  (C.6),  (C.7)  and  (C.  1 1),  whichever  is 
applicable.  In  Case  la,  for  example,  where  layer  1  properties  are  being  modified  by  adding  water 
from  layer  0,  this  means 


*  _  t\(P2  ~P\)  +  T0(pi  -px) 

1 1  — - ~ - 1 

Pi  ~  Pi 

e  _  Si{p2 -p\)  +50(pi  -px) 

Pi  -Pi 

The  new  values  f\ ,  S\  would  combine  to  produce  the  desired  specific  volume  oq  in  layer  1  if  the 
equation  of  state  were  linear,  i.e.,  if 

were  exactly  satisfied.  This  not  being  the  case,  the  specific  volume  implied  by  f\,S\  will  dif¬ 
fer  slightly  from  the  prescribed  value  ofi.  In  cases  where  ofi  equals  the  layer  reference  value  oq,  this 
difference  prevents  the  scheme  from  exactly  restoring  isopycnic  conditions  in  one  instant.  This 
disadvantage  is  likely  to  be  inconsequential,  however,  because  repetitive  use  of  the  coordinate 
restoration  scheme  presented  here  will  quickly  reduce  the  residual  difference  to  zero. 

Hybrid  models  carrying  the  pair  a,  S  as  prognostic  variables  while  diagnosing  T  from  the 
equation  of  state  do  not  suffer  from  this  particular  complication,  because  a  in  such  models  is  not 
affected  by  the  process  of  redistributing  S  among  layers.  Note,  however,  that  the  curvature  of 
isopycnals  in  T,  S  space  is  such  that  a,  S-conserving  models  will  see  a  small  but  systematic  rise  in 
/  T  dp  as  a  result  of  vertical  remapping. 

Division  of  a  coordinate  layer  into  two  sublayers  of  differing  densities  and  extruding  one  for  the 
purpose  of  altering  the  layer’s  original  density  requires  special  consideration.  The  strategy  pres¬ 
ently  adopted  in  Case  2a  is  to  unmix  T  and  S  in  a  direction  perpendicular  to  the  isopycnal  oq  in 
T /S  space.  Details  are  as  follows. 

Define  nondimensional  temperature  and  salinity  as  T  =  T/Tscl  and  S'  =  S/SscX.  Setting  the  left- 
hand  side  in  the  expression 


,  0a  ,  ,  0a  ,  , 

da  =  erdr+0ydx 


to  zero  reveals  that  the  vector  (dT',  dS")  pointing  in  the  direction  tangential  to  an  isopycnal  in  T /S 
space  must  be  perpendicular  to  the  vector  (0a/0T', 0a/05').  Hence,  the  direction  normal  to  the 
isopycnal  is  perpendicular  to  the  vector  (0a /dS',  —0a /dT').  This  can  be  stated  in  the  form 
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dS'  _  d  V 
0a/05"  “  0a/0T' ' 

Combining  the  last  two  expressions,  we  find  that  the  increments  (dr',dS')  changing  the  specific 
volume  a  to  a  +  da  in  a  direction  perpendicular  to  the  isopycnal  a  are  given  by 


dr'  =  da 


8a/er 

(0a/0T')2  +  (da/dS7)2  ’ 


(C.13) 


dS'  =  da 


0a/0<S" 

(0a/0r')2  +  (0a/05')2 


(C.  14) 


Defining  T,S  in  the  two  sublayers  in  Case  2a  is  now  simply  a  matter  of  defining  da  in  (C.13)  and 
(C.  14).  To  set  the  T/S  properties  in  the  lower  sublayer,  we  set  da  =  6q  —  ai  while  in  the  upper 
sublayer  we  set  da  =  ao  —  ai. 

Note  that  “perpendicularity”  in  T/S  space  depends  on  the  scaling  values  TscUSsc j.  The  relative 
contribution  of  T  and  S  unmixing  to  the  creation  of  the  required  density  contrast  between  the 
sublayers  can  be  controlled  by  these  scaling  values.  For  example,  as  Tsd/Ssd  — >  oo,  density  contrast 
is  achieved  entirely  by  T  unmixing. 

With  the  exception  of  Case  2a,  the  density  restoration  method  outlined  above  may  be  viewed  as 
an  upstream  or  donor  cell  advection  scheme  whose  diffusive  character  is  a  possible  concern. 
Schemes  that  modify  the  density  of  both  the  receiving  and  the  donating  layer  (an  attribute  found 
in  the  less  diffusive  centered-differencing  schemes)  offer  a  way  to  reduce  this  numerical  diffusion. 
The  density  restoration  scheme  sketched  below  attempts  to  mimic  the  diffusion  properties  of 
centered  differencing  in  a  manner  similar  to  the  approach  taken  in  Case  2a.  Specifically,  this 
scheme  divides  layer  1  into  two  sublayers,  one  of  which  has  the  desired  target  density  ofi.  The 
second  sublayer  is  then  expelled  into  either  the  layer  above  or  below,  depending  on  whether  the 
second  sublayer  is  lighter  or  denser  than  the  first  one. 

Of  utmost  importance  in  the  design  of  such  an  “unmixing”  scheme  is  that  it  should  not  create 
artificial  T,S  extrema  in  the  water  column.  To  guarantee  this,  we  require  the  unmixing  process  to 
be  confined  to  a  “bounding  box”  in  T /S  space  whose  four  sides  are  defined  in  terms  of  T,S  values 
in  layers  0  and  2.  Since  the  purpose  of  unmixing  is  to  create  sublayers  of  prescribed  density,  one  or 
two  corners  of  the  bounding  box  may  in  practice  be  truncated  by  isopycnals  representing  the 
target  densities  of  the  two  sublayers. 

The  sublayer  generation  routine  presently  in  use  treats  the  target  density  and  T,  S  con¬ 
straints  as  independent  and  requires  that  both  be  satisfied.  The  scheme  tests  two  “candidate” 
unmixing  directions  as  to  their  ability  to  create  sublayers  of  the  desired  densities.  Both  un¬ 
mixing  directions  are  defined  as  lines  that  connect  the  T,S  point  representing  the  original  layer 
1  in  T/S  space  with  either  the  upper  left  or  lower  right  corner  of  the  bounding  box.  In  either 
case,  the  unmixing  line  terminates  on  both  ends  on  either  the  target  isopycnal  or  the  edge  of 
the  bounding  box,  whichever  is  encountered  first.  The  unmixing  direction  ultimately  chosen  is 
the  one  that  results  in  the  bigger  density  contrast.  Knowing  the  end  points  of  the  unmixing 
line,  we  can  compute  the  thickness  of  the  two  sublayers  by  requiring  T,S  conservation  during 
unmixing. 

Examples  of  unmixing  scenarios  are  shown  in  Fig.  12.  Abscissa  and  ordinate  in  the  three  plots 
are  salinity  and  temperature,  respectively.  The  bounding  box  is  dashed  while  the  slanting  curves 


R.  Bleck  /  Ocean  Modelling  37  (2002  )  55-88 


85 


Fig.  12.  Three  examples  of  layer  unmixing  constrained  by  T/S  bounding  boxes  (dashed)  and  bounding  isopycnals 
(gently  curved  lines).  Filled  and  hollow  diamonds  mark  the  T /S  values  before  and  after  unmixing,  respectively. 


represent  isopycnals.  The  black  diamond  represents  the  layer  to  be  unmixed,  while  the  hollow 
diamonds  represent  the  two  resulting  sublayers. 

Aside  from  reducing  numerical  diffusion,  this  unmixing  scheme  removes  an  undesirable  bias 
from  the  coordinate  restoration  logic,  namely,  that  layers  subject  to  density  restoration  invariably 
grow  at  the  expense  of  their  neighbors.  The  present  unmixing  scheme  has  the  opposite  effect,  i.e.,  it 
causes  layers  whose  density  is  being  manipulated  to  shrink.  We  use  the  cushion  function  described 
earlier  to  prevent  a  layer  from  getting  too  thin  in  the  process.  Our  present  strategy  is  to  attempt 
partial  coordinate  restoration  through  unmixing  before  resorting  to  the  donor  cell  method  de¬ 
scribed  earlier.  The  exact  meaning  of  the  word  “partial”  remains  to  be  determined. 

Since  the  grid  generator  consumes  little  time  compared  to  other  parts  of  HYCOM,  timing  tests 
fail  to  show  a  systematic  impact  of  the  unmixing  scheme  on  model  performance. 

In  the  present  grid  generator  implementation,  the  parameter  <5  appearing  in  (C.4)  and  (C.5)  is 
set  to  10  m.  Except  in  very  shallow  water,  it  is  advisable  to  make  <5  a  constant.  This  will  prevent 
the  potential-vorticity  conserving  advection  and  Coriolis  terms  in  the  momentum  equations, 
which  are  written  in  a  shallow-water  form  inherited  from  MICOM,  from  generating  spurious 
dynamic  stretching  effects.  Keeping  coordinate  surfaces  horizontal  in  geographic  regions  where 
layers  are  unable  to  retain  their  isopycnal  character,  i.e.,  where  coordinate  surfaces  are  not  ma¬ 
terial  under  adiabatic  flow  conditions,  assures  that  layer  thickness  drops  out  of  the  momentum 
advection  and  Coriolis  terms,  allowing  them  to  act  as  they  would  in  Cartesian  coordinate  models. 

In  shallow  coastal  areas,  on  the  other  hand,  it  is  appropriate  to  scale  <5  by  the  total  water  depth. 
The  argument  here  is  that  Taylor  columns  in  the  coastal  zone  are  likely  to  span  the  whole  water 
depth,  implying  that  undulations  in  the  depth  of  er  surfaces  (cr  =  depth  scaled  by  bottom  depth) 
accurately  portray  the  stretching  effect  due  to  sea  floor  undulations. 


Appendix  D.  Transformation  of  nonisopycnic  model  output  to  isopycnic  coordinates 

Let  p[p)  be  a  piecewise  constant  (i.e.,  stairstep)  density  profile  where  the  height  A p  of  each 
stairstep  (“riser”)  represents  the  thickness  of  an  isopycnic  layer  while  the  width  A p  of  each 
stairstep  represents  the  density  increment  between  consecutive  isopycnic  layers.  Both  step  widths 
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and  riser  heights  are  arbitrary.  Our  task  is  to  transform  this  profile  into  another  stairstep  profile 
differing  from  the  original  one  in  that  the  location  of  the  risers  on  the  density  axis  is  prescribed. 
The  transformation  should  preserve  the  mean  density  in  the  column,  i.e.,  the  integral  f  p  dp  taken 
from  the  sea  surface  to  the  bottom. 

Let  pk(k  =  1 pk+ j  >  pk)  mark  the  points  on  the  density  axis  where  we  want  the  new  risers 
to  be  placed.  We  require  that  the  pk  span  the  density  range  of  the  input  profile,  and  that  the  input 
profile  be  monotonic.  Denoting  the  height  of  the  upper  and  lower  interface  bounding  the  /cth  layer 
by  Pk-i/i  and  pk+\/2,  respectively,  the  condition  we  wish  to  satisfy  can  then  be  stated  as 

"  fP"+ 1/2 

}^Pk(Pk+l/2  -Pk-X/l)  =  /  pdp.  (D.l) 

k=l  Jp\/2 


The  interface  depths  are  the  unknowns  in  the  problem. 

Integration  of  (D.l)  by  parts  -  on  the  left  this  amounts  to  reordering  the  terms  under  the 
summation  sign  -  allows  us  to  rewrite  (D.l)  as 

»-!'  rP(Pn+ l  ft) 

P nPn+\/2  -  PlPl/2  ~  ^Pk+l^Pk+l  ~  Pk)  =  [pp]^  ~  /  pdp. 

k=  1  '  J  piPlft) 

By  virtue  of  p{p\/2)  ^  P\  and  p(pn+ 1/2)  ^  pn,  this  expression  reduces  to 

(  ppn 

Z2Pk+\/i{pk+x  ~  Pk)  =  /  P^P-  (D-2) 

k=  1  J  Pi 


Our  strategy  is  to  satisfy  (D.2)  by  breaking  the  integral  into  pieces  taken  over  intervals  (pk,pk+  x) 
and  conserving  each  integral  individually.  This  immediately  leads  to 

1  rPk+i 

Pk+ 1/2= - - —  /  pdp  (k=  1). 

Pk+ 1  Pk  J  pk 

The  procedure  for  transforming  the  density  profile  is  illustrated  in  Fig.  13. 


Pk-2  Pk-1  Pk  Pk+1 


P 


Fig.  13.  Illustration  of  vertical  transform  procedure.  Thin/thick  stairstep  lines  represent  input/output  density  profiles. 
Requiring  the  hatched  areas  to  be  of  equal  size  yields  location  of  pk-i/i-  See  text  for  additional  details. 
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Knowing  the  new  interface  depths,  we  can  now  proceed  to  transform  other  model  output 
variables.  Care  must  be  exercised  in  transforming  mass  fluxes  (that  is,  the  product  of  velocity  and 
layer  thickness),  because  their  vertical  sum  needs  to  be  preserved  if  they  are  to  be  used  for 
quantitative  diagnostic  work  in  density  space.  Note  that  preserving  the  vertical  sum  of  mass  fluxes 
is  equivalent  to  preserving  the  vertical  integral  of  the  associated  velocity  field.  Hence,  if  constraints 
patterned  after  (D.l)  are  used  in  transforming  mass  fluxes,  the  variable  replacing  p  in  (D.l)  must 
be  velocity,  not  mass  flux. 

In  layer  models  like  MICOM  and  HYCOM,  where  mass  flux  computation  is  a  multi-step 
process  involving  a  spatially  varying  mix  of  diffusive  and  antidiffusive  fluxes,  the  first  step  in  the 
flux  transformation  therefore  is  generation  of  an  auxiliary  velocity  field  representing  mass  flux 
divided  by  layer  thickness.  It  is  this  velocity  field  that,  after  being  transformed  in  a  manner 
preserving  its  column  integral,  forms  the  basis  for  mass  fluxes  in  pk  space. 
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